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ABSTRACT 

We analyze the frequency evolution of millisecond oscillations observed during type I X-ray bursts 
with the Rossi X-ray Timing Explorer in order to establish the stability of the mechanism underlying 
the oscillations. Our sample contains 68 pulse trains detected in a search of 159 bursts from 8 accreting 
neutron stars. As a first step, we confirm that the oscillations usually drift upward in frequency by about 
1% toward an apparent saturation frequency. Previously noted anomalies, such as drifts toward lower 
frequencies as the oscillations disappear ("spin-down" episodes) and instances of two signals present 
simultaneously at frequencies separated by a few Hz, occur in 5% of oscillations. Having verified the 
generally accepted description of burst oscillations, we proceed to study the coherence of the oscillations 
during individual bursts, and the dispersion in the asymptotic frequencies in bursts observed over five 
years. On short time scales, we find that 30% of the oscillation trains do not appear to evolve smoothly 
in phase. This suggests either that two signals are present simultaneously with a frequency difference 
too small to resolve (< 1 Hz), that the frequency evolution is discontinuous, or that discrete phase jumps 
occur. On time scales of years, the maximum frequencies of the oscillations exhibit fractional dispersions 
of A^max/ (fmax) £4x 10~ 3 . In the case of 4U 1636—536, this dispersion is uncorrelated with the known 
orbital phase, which indicates that a mechanism besides orbital Doppler shifts prevents the oscillations 
from appearing perfectly stable. In the course of this analysis, we also search for connections between 
the properties of the oscillations and the underlying bursts. We find that the magnitudes of the observed 
frequency drifts are largest when the oscillations are first observed at the start of the burst, which suggests 
that their evolution begins when the burst is ignited. We also find that radius expansion appears to 
temporarily interrupt the oscillation trains. We interpret these results under the assumption that the 
oscillations originate from anisotropies in the emission from the surfaces of these rotating neutron stars. 

Subject headings: stars: neutron — X-rays: bursts — X-rays: stars 



1. INTRODUCTION 

Millisecond oscillations have been observed during ther- 
monuclear X-ray bursts from nine neutron stars in low 
mass X-ray binaries (LMXBs; see Strohmayer 2001, for a 
review) . The bursts are triggered by the unstable nuclear 
burning of accreted material on the neutron star's sur- 
face (see Lewin, van Paradijs & Taam 1993, for a review). 
It therefore has long been expected that anisotropies in 
the burning should produce pulsations at the stellar spin 
frequency (e.g Bildsten 1995; Strohmayer et al. 1996). 
Indeed, the general characteristics of the observed os- 
cillations suggest that they originate from a brightness 
anisotropy on the stellar surface. First, the oscillations are 
highly coherent towards the end of the burst (Strohmayer 
et al. 1996; Strohmayer & Markwardt 1999). Second, the 
oscillations are frequently observed in the rise of bursts, 
when there is spectral evidence for a growing burning re- 
gion (Strohmayer, Zhang, & Swank 1997). Third, the 
frequencies of these oscillations are remarkably similar in 



bursts separated by several years (Strohmayer et al. 1998; 
Muno et al. 2000), which suggests that they originate from 
a stable clock. 

If the burst oscillation frequency, t'burst, is indeed the 
spin frequency of the neutron star, the oscillations are key 
to understanding several aspects of these systems. Pairs 
of kilohertz quasi-periodic oscillations (kHz QPOs) are ob- 
served in the persistent emission from many neutron star 
LMXBs. Their frequency difference AfkHz remains com- 
parable to 0.5 or 1 x i^burst, even though the frequencies of 
the kHz QPOs vary by more than a factor of 2 (see van 
der Klis 2000, for a review). It has been suggested that 
these QPOs result from signals at a Keplerian frequency 
in the disk and a beat between that frequency and the spin 
of the neutron star (Strohmayer et al. 1996; Miller, Lamb, 
& Psaltis 1998). 2 If this is the case, we can infer the ro- 
tation rate of the neutron stars in more than 20 systems 
(van der Klis 2000). The inferred spin frequencies cluster 
around 300 Hz. This confirms the suggestion that these 
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2 Alternative models for the kHz QPOs cannot account naturally for the similarity between i-'burst an< i ^kHz ( e -g Psaltis 2001). 
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neutron star LMXBs are the progenitors of recycled mil- 
lisecond radio pulsars (Alpar et al. 1982; Radhakrishnan 
& Srinivasan 1982). However, it also implies that some 
mechanism limits the frequencies to which these stars can 
be spun-up by accretion, such as a correlation between 
the surface magnetic field strength and X-ray luminosity 
(White & Zhang 1997) or gravitational radiation from a 
quadrupole in the mass distribution of the star (Bildstcn 
1998). 

A closer examination of these burst oscillations reveals 
an even more complex picture (see also Psaltis 2001, for 
a detailed discussion). First, the coherence of the oscilla- 
tions is cast into some doubt by the detection of a sudden 
0.25 cycle phase shift in one oscillation out of seven stud- 
ied by Strohmayer (2001). Second, oscillations are not 
only observed in the rise of thermonuclear bursts, but also 
far into the decay (e.g. Smith, Morgan, & Bradt 1997), 
even though by this time the fuel over the entire surface of 
the neutron star is expected to have been consumed by the 
burning front leaving no observable brightness anisotropics 
(e.g. Fryxell & Woosley 1982). Third, while the spin of 
the neutron star must be constant during the course of 
a burst, the frequency of the oscillations is usually ob- 
served to increase and saturate at some maximum fre- 
quency (Strohmayer et al. 1997). It is this "asymptotic" 
frequency that appears stable in bursts separated by sev- 
eral years. The frequency drift is thought to originate from 
one of several mechanisms that might carry anisotropies in 
a retrograde motion about the neutron star, resulting in 
an observed frequency slightly lower than that of the star's 
spin (Strohmayer et al. 1997; Cumming & Bildsten 2000; 
Spitkovsky, Levin, and Ushomirsky 2002; Heyl 2002). Fi- 
nally, in a handful of bursts there are episodes of drifts to 
lower frequency ("spin-down") in the tails of the bursts 
(Strohmayer 1999; Miller 2000; Muno et al. 2000) and 
strong detections of signals at second frequencies within 
a few Hz of the primary signal (Miller 2000; Galloway et 
al. 2001). 

If anything is clear, it is that more observational clues 
are required to understand burst oscillations. In this pa- 
per, we address the how the frequency of these oscillations 
evolves as a function of time, in order to establish the co- 
herence of the underlying mechanism. In particular, we 
employ standard pulse phase analysis techniques adapted 
from pulsar timing studies (Manchester & Taylor 1977). 

2. OBSERVATIONS AND DATA ANALYSIS 

Our analysis used observations with the Proportional 
Counter Array (PCA; Jahoda et al. 1996) on the Rossi X- 
Ray Timing Explorer (RXTE). PCA consists of five iden- 
tical gas-filled proportional counter units with a total ef- 
fective area of 6000 cm 2 and sensitivity in the 2.5-60 keV 
range. The detector is capable of recording photons with 
microsecond time resolution and 256-channel energy reso- 
lution. The data were recorded in a wide variety of data 
modes with different time and energy resolutions, depend- 
ing upon the details of the original proposed programs and 
the available telemetry bandwidth. For all of the analysis 
presented here, we converted the TT photon arrival times 
at the spacecraft to Barycentric Dynamical Time (TDB) 

3 Oscillations in bursts associated with MXB 1743—29 have been observed during observations of the bursting pulsar GRO J1744— 28. A search 
for these bursts was not part of the analysis presented here. 



Fig. 1. — A dynamic power spectrum illustrating the typical fre- 
quency evolution of a burst oscillation. Contours of power as a 
function of frequency and time were generated from power spectra 
of 2 s intervals computed every 0.25 s. A Welch function was used to 
taper the data to reduce sidebands in the power spectrum due to its 
finite length (Press et al. 1992). The contour levels are at powers 
of 0.02 in single-trial probability starting at a chance occurrence of 
0.02. The PCA count rate is plotted referenced to the right axis. 

at the solar system barycenter, using the Jet Propulsion 
Laboratory DE-200 solar system ephemeris (Standish et 
al. 1992).' 

We searched the entire RXTE public data archive for X- 
ray bursts from 8 neutron stars 3 that are known to exhibit 
burst oscillations (see Table 1 and Muno et al. 2001). As 
of September 2001, we have identified a total of 159 X-ray 
bursts from these 8 sources. Each of these bursts was then 
searched for millisecond oscillations as described below. 

For our analysis, we used data containing a single en- 
ergy channel (2.5-60 keV) and 2 -13 s (122 /xs) resolution. 
We note that the use of more restricted energy bands can 
sometimes allow oscillation trains to be detected for a sec- 
ond or two longer, since the oscillation amplitudes are of- 
ten stronger at higher energies (e.g., Muno et al. 2000; 
Giles et al. 2002). However, in many cases, the data modes 
with higher energy resolution contained data gaps due 
to the limited size of the satellite memory buffers, mak- 
ing it difficult to produce a coherent phase model span- 
ning the entire oscillation train (see Section 2.2). Data 
modes specifically designed to start recording during a 
burst ("burst catcher" modes) were analyzed in combi- 
nation with the other modes in order to fill data gaps, but 
the burst catcher modes with high time resolution only 
recorded in a single 2.5-60 keV energy channel. 

2.1. Fourier Timing 

We produced Fourier power spectra of 1 s intervals of 
data for the first 16 s of each burst, and searched for signals 
within a frequency range of ±5 Hz from the frequencies 
listed in Table 1. Any oscillation with a single trial proba- 
bility of less than 2 x 10 -5 that it is due to noise (a power 
of P > 21.6 according to the normalization of Leahy et al. 
1983) was considered a detection. This corresponds to a 
1% probability that a noise signal will be that strong given 
a search of an entire burst, so that 1-2 noise signals may 
appear significant from the search of our entire sample of 
bursts. We detected 68 oscillation trains in this manner 
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(Tabic 1). We also recorded the frequencies of each signal 
detected. The frequencies in Table 1 were determined re- 
cursively, by computing the median value (rounded to the 
nearest Hz) from the signals detected until that value did 
not change. The range of detected frequencies is reported 
in Section 3.1. 

To qualitatively examine the frequency evolution, we 
computed power spectra of 2 s intervals of data every 0.25 
s, and plotted contours of power as a function of both 
time and frequency. We tapered the data with a Welch 
function (Press et al. 1992) before computing each power 
spectrum, in order to suppress the power in sidebands due 
to the window function. A typical example of such a dy- 
namic power spectrum is illustrated in Figure 1. We use 
dynamic power spectra to characterize the oscillations in 
Sections 3.2-3.3. 

2.2. Phase Connection 

In order to examine in detail whether the oscillations 
are coherent and drift to stable asymptotic frequencies, 
we developed frequency models using the phase connec- 
tion method described in Muno et al. (2000). We modeled 
all oscillations that lasted continuously for more than 2 s 
using the 2~ 13 s data (Table 1). A summary of the proce- 
dure is as follows. We fold the data in 0.25-0.5 s intervals 
according to a trial frequency model (which is the deriva- 
tive of a phase model), and compute a Fourier transform 
of the resulting profile. If the power P in the fundamental 
frequency of the pulsation 4 has less than a 10% chance of 
occurring randomly due to noise (P > 4.6 in the normal- 
ization of Leahy et al. 1983), we measure its phase by a 
linear least squares fit of a sinusoid to the profile. The 
measured phase residuals </> bs(ii) are then related to the 
predicted phases mo dei(ii) by 

A<j>(ti) = 4>obs{ti) - ^model(ij)- (1) 

In the above convention, an increasing phase residual indi- 
cates that the instantaneous trial frequency (d(f> mo <i e \/ dti) 
is too low. 

The phase residuals A<f>(U) are then modeled via a least 
squares minimization, using a x 2 statistic calculated from 
the phase residuals: 



This procedure can be iterated using the new frequency 
model, if it is composed of linear functions. 

The advantages of this method are that (i) the form of 
the observed phase residuals naturally suggests the best 
model for the frequency evolution and provides a quanti- 
tative test of how well the trial phase evolution matches 
the data, (ii) there exist standard techniques to calculate 
the uncertainties on the model from the least squares fit to 
the phases, and (Hi) the method utilizes the phase informa- 
tion that the power spectrum discards, which allows much 
finer frequency resolution on short time scales. In con- 
trast, the Z 2 method of Strohmayer & Markwardt (1999) 
distinguishes between models based upon the amount of 
power in the folded profile, but provides little information 
about how the frequency of an oscillation deviates from 
the model assumed. The Z 2 statistic is still very useful 
for examining the amplitude and phase evolution of the 
oscillations on short (less than 0.1 s) time scales after the 
frequency evolution is modeled. We apply this technique 
and use the resulting models in Sections 3.4-3.7. 

2.3. Energy Spectra 

Finally, we produced energy spectra in 0.25 s inter- 
vals from each burst using available combinations of data 
modes that provide at least 32 energy channels. The de- 
tector response was estimated using PCARSP in FTOOLS 
version 5.1 (see http://heasarc.gsfc.nasa.gov/lheasoft/). 
We subtracted spectra from 15 s of emission prior to the 
burst to account for background, and fit each spectrum 
between 2.5-20 keV with a model consisting of a black- 
body multiplied by a constant interstellar absorption fac- 
tor. The absorbing column for each burst was taken to be 
the mean value from fits to each 0.25 s interval assuming a 
variable absorption. The model provides a color tempera- 
ture (T co i) and a normalization proportional to the square 
of the apparent radius of the burst emission surface (R co i), 
and allows us to estimate the bolometric flux as a function 
of time. We compare these spectra to the properties of the 
oscillations in Sections 3.3 and 3.6. 

3. RESULTS 

3.1. Frequency Range of Oscillations 

In order to establish the range of frequencies over which 
the oscillations are observed, and in particular whether 
there are strict upper limits to their frequencies, we pro- 
duced a histogram of the distribution of the frequencies 
at which signals were detected in non-overlapping 1 s in- 
tervals during the bursts from each source (Figure 2). 
For the majority of sources, the histograms are narrowly 
distributed with a width of 2-4 Hz, significantly smaller 
than the 10 Hz search window. This corresponds to fre- 
quency drifts ranging from 0.4% (KS 1731-260) to 1.2% 
(4U 1702-429). 

The distribution from 4U 1636—536 appears to be 
double-peaked, possibly owing to two effects. First, os- 
cillations are usually observed from this source only in 
the rise and tail, and the frequency gap between the two 
peaks could result from the fact that we do not observe 
the intermediate frequencies as the oscillations evolve (see 

4 We have confirmed that all of the profiles were consistent with sinusoidal signals (Muno et al. 2000, 2001), justifying the use of only the 
power at the fundamental frequency. 
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Here, o~A4> is the uncertainty on the phase measured from 
the folded profile, which can easily be determined from the 
diagonal values of the covariance matrix in the linear least 
squares fit (Press et al. 1992). We have confirmed that 
these uncertainties are accurate via Monte Carlo simula- 
tions, in which we (i) perturbed the counts in each bin of 
the profile within the range expected from Poisson count- 
ing noise and (ii) recorded the uncertainty as the standard 
deviation of phase measurements of 100 perturbed profiles. 
We also verified the values for cta<£ by visual inspection. 

The phase connection technique provides a correction 
d(A(j))/dt to the original frequency model f mo dci, so that 
the best-fit frequency is 
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Fig. 2. — Histograms of the number of oscillations detected at each frequency within a 10 Hz range of the mean indicated during the first 
16 s of each burst, separated by source. Detections are defined as having a single-trial probability of less than 2 X 10~ 5 that they are due to 
noise. The histograms arc fairly narrowly distributed with an absolute width of about 4 Hz. Note that the signal at 571 Hz in MXB 1659—298 
has a chance probability of 3xl0 -4 given that we searched 159 bursts from all sources to produce the histograms. 



Section 3.2). However, two oscillations are present simul- 
taneously in several bursts from this source (see 
Section 3.3 and Miller 2000), which could indicate there 
are indeed two fundamental frequencies. 

One detected oscillation is unusual given the distribu- 
tions in Figure 2, and has already been noted by Wij- 
nands, Strohmayer, & Franco (2000). It is detected in 
MXB 1659-298 about 4 Hz above the next highest fre- 
quency, with a probability of 4 x 10~ 4 that it is due to 
noise given all of the trial frequencies searched. The oscil- 
lation occurs at 571 Hz in a single interval about 15 s into 
the tail of a burst, while the oscillation detected earlier 
in the burst seemed to disappear about 5 s into the burst 
after reaching an asymptotic frequency of 567 Hz. If this 

signal is a burst oscillation, it calls into question whether 
we are observing frequency saturation at all. 

3.2. Instances of Unusual Frequency Evolution 

We also searched the dynamic power spectra to con- 
firm that the oscillations in our sample drift upward in 
frequency during the course of the burst, saturating at an 
approximately constant frequency late in the tail. This 
trend occurs in most bursts, and is illustrated in Figure 1. 
However, there are a few exceptions that have been noted 
in the literature, such as (i) clear evidence for a frequency 
decrease several seconds after the burst is first detected 



(spin-down) and (ii) signals present at two frequencies si- 
multaneously. Spin-down episodes and simultaneous sig- 
nals are not common. 

We find only 3 clear examples of spin-down episodes dur- 
ing bursts, which we display in Figure 3. Two of these have 
been reported previously, from 4U 1636—536 (Strohmayer 
1999) and KS 1731-260 (Muno et al. 2000), and we find 
one additional example from 4U 1728—34. In all cases, 
the frequency decrease is on order 1 Hz. There is no clear 
pattern as to when in the burst the spin-down occurs. In 
KS 1731—260, the event occurs in the peak of the burst, 
while in 4U 1728-34 and 4U 1636-536 the spin-down oc- 
curs as the burst decays. 

We also find two clear examples of bursts in which two 
signals are detected simultaneously at frequencies sepa- 
rated by about 1 Hz. The most significant example has al- 
ready been reported by Miller (2000) from 4U 1636-536, 
and we find a second example in 4U 1702—429. Both 
are displayed in Figure 4. The second oscillation in 
4U 1636—536 occurs between 2-3 s into the burst, and 
has a chance probability < 3 x 10~ 10 of occurring ran- 
domly due to noise (see also Miller 2000). The second 
signal from 4U 1702—429 occurs 11 s into the burst, and 
is only marginally significant, with a single-trial chance 
probability of 5 x 10~ 5 . Less significant examples, with 
chance probabilities of ~ 0.01, have also been reported in 
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Fig. 4. — Same as Figure 1, for oscillations that exhibited two 
simultaneous signals separated in frequency by ~ 1 Hz. The sec- 
ondary signal occurs at 2—3 s in panel a, with a chance probability 
that it is due to noise of 3 X 10~ 10 . The second signal at 11 s in 
panel b has a chance probability of 5 X 10~ 5 . 

be examined, because there were gaps in the high time 
resolution data during the peak. Since it also has been no- 
ticed previously that oscillations in KS 1731—260 (Smith, 
Morgan, & Bradt 1997) and MXB 1743-29 (Strohmayer 
et al. 1997) appear only after radius expansion ends, we 
conclude that radius expansion either obscures or inter- 
rupts the oscillations. 



Fig. 3. — Same as Figure 1, for the 3 oscillations that exhibited 
spin-down episodes out of our sample of 68 pulsations. The fre- 
quency drift is on order 1 Hz in each case. There appears to be no 
trend as to when in the burst spin-down occurs. 

two more bursts from 4U 1636—536 by Miller (2000), and 
in one burst from 4U 1916-053 by Galloway et al. (2001). 
In all cases, the secondary frequencies are 1-3 Hz lower 
than the main signal. 

3.3. How Radius Expansion Affects Oscillations 

In many bursts, photosphcric radius expansion is evi- 
dent at the start of the burst, during which _R co i increases 
and T co i decreases such that the bolometric flux remains 
constant, presumably at the Eddington limit (see Lewin et 
al. 1993). We plot examples of bursts with and without ra- 
dius expansion in Figure 5. We examined whether radius 
expansion affects the times in the burst when the oscil- 
lations are observed. We find that all 16 of the bursts in 
which oscillation trains are observed continuously through- 
out the rises, peaks, and tails fail to exhibit radius ex- 
pansion (e.g. Figure 5a). Of the 51 bursts where the 
oscillation train is not observed during the peak, 40 ex- 
hibit radius expansion (e.g. Figure 5b). One burst from 
4U 1608-522 (1998 March 27 at 14:08:30 TDB) could not 



3.4. Mathematical Form of the Frequency (Phase) 
Evolution 

We modeled the frequency evolution of the oscillations 
using the phase connection technique described in Sec- 
tion 2.2. We tried both polynomial models (up to fifth or- 
der) and a saturating exponential such that the frequency 
evolution is of the form used by Strohmayer & Markwardt 
(1999): 

v{t) = v Q + Ai/exp(-t/r). (4) 

Here, z^o is the asymptotic frequency of the oscillation, Av 
is the frequency drift referenced to the start of the burst, 
and t is the time scale for the frequency drift. When mod- 
eling the data with polynomials, we iterated the fitting 
procedure until no improvement in the fits to the phase 
residuals was found. The fits using the exponential model 
were executed on phase residuals found assuming a con- 
stant trial frequency, as the corrections must add linearly 
to the trial phase model in order to use the iterative pro- 
cess. Table 2 lists the number of oscillations best fit with 
each model, as determined from the value of reduced \ 2 . If 
the reduced x 2 was greater than one, and if both the num- 
ber of degrees of freedom and the value of % 2 decreased in 
comparing one model to another, an F-test (Bevington 
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Fig. 5. — Dynamic power spectra compared to parameters from blackbody fits to the energy spectra. Top Panels: Same as Figure 1. 
Center Panels: The color radius derived from the spectral fits, as a function of time. Bottom Panels: The color temperature of the burst. 
Uncertainties arc l-cr. No radius expansion occurred during the burst on the left, and the oscillations were observed continuously for the first 
7 s of the burst. The radius expansion on the right appears to interrupt the oscillations. Note that the power spectra are of 2 s intervals 
spaced every 0.25 s, while the energy spectra are computed in non-overlapping 0.25 s intervals. 
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Fig. 6. — The frequency and phase evolution of burst oscillations from (a) 4U 1636—536 and (b) 4U 1702—429. The figures demonstrate 
how the phase connection method allows one to distinguish between models for the frequency evolution. Top panel: Same as Figure 1. The 
solid line at the center of the contours is the best-fit frequency evolution. Center and bottom panels: The phase residuals after folding a the 
data about the model for the frequency evolution indicated. The bottom panel is the best-fit, as indicated by the lower value of xt ■ 



1969) was used to ensure that the decrease in \ 2 had less 
than a 5% chance of being random. 

Figure 6a illustrates the importance of considering the 
phases in modeling the frequency evolution. Judging from 
power spectra alone, the frequency of this oscillation from 
4U 1636—536 appears roughly constant (top panel). How- 
ever, when one examines the phase residuals from a con- 
stant frequency model it is clear that additional evolution 
has occurred (center panel). An acceptable fit to the fre- 
quency evolution only can be achieved using a fifth or- 
der polynomial (bottom panel). Although some of this fre- 
quency evolution can be inferred from the power spectra 
alone (see also Miller 2000), by tracking the phase we are 
able to achieve high frequency resolution and to determine 
how coherent the oscillations are. 

We modeled 59 oscillations from 6 different sources (Ta- 
ble 1). The results are summarized in Table 2. Thirty- 
seven oscillations required either an exponential model or 
a polynomial of third degree or higher. Exponential mod- 
els are only favored over polynomials in only 6 out of 37 
cases. The exponential models are formally inconsistent 
(greater than 99% confidence) in 22 of the remaining 31 
cases, while all such frequency models are inconsistent with 
the data in 12 out of 37 cases. 

The polynomial models generally reproduce the fre- 
quency evolution better than exponential models because 
there is slower frequency evolution at the start of the oscil- 
lation train, inflection in the evolution during the train, or 



a decrease in the oscillation frequency as the trains disap- 
pear (see also Miller 2000). Figure 6 exhibits all of these 
variations in the frequency evolution. In addition to the 
upward drifts in frequency that saturate after a few sec- 
onds, the frequencies also tend to wander by about 0.1 Hz 
(see, e.g., the evolution in the top panel of Figure 6a). 

We see no evidence for a rapid decrease in frequency at 
the start of the oscillation that is comparable in magni- 
tude to the later upward drift in frequency Such a fre- 
quency decrease could occur as the burning front spreads 
around the surface of the neutron star (e.g., Strohmayer 
et al. 1997; Spitkovsky et al. 2002), but may be too rapid 
to observe using our analyses. 

3.5. Short-term Stability of the Burst Oscillations 

In general, we find many more large values of \ 2 than 
would be expected if the phase evolution were described by 
smooth, low-order polynomial and exponential functions. 
In particular, 30% of the oscillation trains are inconsis- 
tent with smooth models with 99% confidence. This could 
indicate that the underlying mechanism is not perfectly 
stable. We tried higher-order polynomials in instances 
where our best-fit model is statistically unacceptable, but 
we find that we achieve only marginal improvements in 
X 2 ■ Higher-order polynomials can smooth residuals in the 
middle of the oscillation train, but tend to create larger 
residuals at the endpoints. Figure 7 illustrates two exam- 
ples of oscillation trains that we were unable to fit with 
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O 5 10 15 5 10 15 

Time (s) since 05:16:42 (TDB) Time (s) since 23:21:59 (TDB) 



Fig. 7. — The frequency and phase evolution of oscillations from two trains for which the best-fit model was statistically inconsistent with 
the data. The axes are as in Figure 6. 



exponential or polynomial models. Piecewise smooth 
phase functions would provide better fits in some in- 
stances, such as between 2-4 s in Figure 7a, and 3-6 s 
in Figure 7b. However, the resulting discontinuities in the 
derivatives of the phases would imply that the frequencies 
of the oscillations shift in less than 0.25 s. 

Moreover, some of the residuals appear to be discrete 
changes in the phases of the oscillations. One possible ex- 
planation for the phase residuals in the first 1-2 s of Fig- 
ure 7a is that two signals are present simultaneously within 
about 1 Hz, the resolution of the folding technique. Simul- 
taneous signals with only slightly larger frequency separa- 
tions have been seen in other bursts from 4U 1636—536 
(Miller 2000). In Figure 7b, there appear to be discrete 
jumps of about 0.1 cycle in the phase of the oscillation 
from 7-10 s into the burst, similar to that observed by 
Strohmayer (2001) during an oscillation in the rise of a 
burst from 4U 1636—536. We confirmed that the sudden 
changes in the phases of the oscillations arc not due to 
variations in the pulse profile, which is always sinusoidal. 
The presence of phase jumps suggests that the oscillations 
are not strictly coherent. 

3.6. Magnitude and Time Scale of Frequency Drift 

Our models also allowed us to compare the frequency 
evolution to other properties of the bursts. We first mea- 
sured the observed change in frequency Av for the con- 
tinuous portion of each oscillation. We then measured the 
time scale r that represents how long the oscillation took 
to evolve in frequency by (1 — e _1 )Azv = 0.632 x Av from 



the minimum observed frequency (compare Equation 2). 
As measures of the burst time scales, we have fit the de- 
cay in flux of each of these bursts with either one or two 
exponential functions (exp[— t/t^])- We also computed the 
peak bolometric flux and the fluence of each burst, using 
the exponential fit to determine the flux contribution at 
late times. We find that neither IS.u nor r are correlated 
with the above burst properties. 

We have also compared Av and r with the duration of 
the oscillations and the time at which they began. We 
find that r is not correlated with the start time or dura- 
tion of the oscillations. Instead, r is distributed uniformly 
between 0.75-5.5 s. We plot the frequency drift Av as a 
function of these parameters in Figure 8. Oscillations that 
are first observed early in the burst evolve over a larger 
frequency range {top panel), but their frequency drift is 
not correlated with how long the oscillations last {bottom 
panel). This suggests that the evolution is initiated by the 
start of the burst. 

Cummings & Bildstcn (2000) previously noticed that 
the fractional frequency drifts Av/v in sources with ^600 
Hz oscillations are a factor of two smaller than those in 
^300 Hz oscillations. We have indicated the data from 
the faster (^600 Hz) oscillations in Figure 8 with filled 
circles. The continuous portions of the fast oscillations 
tend to start later in the burst (with only two exceptions) 
and to drift by a smaller absolute frequency range. We be- 
lieve that this is due to the fact that the faster oscillations 
tend to occur in bursts with radius expansion (Muno et 
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Start Time of Oscillations (s) 




5 10 15 

Duraiion of Oscillaiion (s) 



Fig. 8. — The total frequency drift in the continuous portion of 
the oscillation trains plotted as a function of the (top panel) the time 
at which the continuous oscillation train began, and (bottom panel) 
the duration of the continuous train. The total frequency drift is 
largest when the oscillations are first observed at the start of the 
bursts, but is not correlated with the duration of the oscillations. 
The solid circles are data from fast (~600 Hz) oscillations, while the 
open circles are from slow (~300 Hz). 

al. 2001), which prevents the oscillations from being ob- 
served earlier in the bursts (see Section 3.1). Indeed, the 
fast oscillators in Figure 2 do not in general have a nar- 
rower distribution of observed frequencies when consider- 
ing the non-continuous portions of the oscillations as well. 
However, we do confirm that the fractional frequency drift 
is a factor of two smaller in the ^600 Hz oscillators than 
in the ~300 Hz ones, but only because the frequencies v 
are larger. 

If the frequencies of the oscillations are associated with 
underlying stable clocks, then the fact that they evolve in 
frequency implies that the phases of the oscillations drift 
by several cycles with respect to those clocks. If we take 
the maximum observed frequency (i/ max ) to represent the 
stable clock (see Section 3.7), then we can define a phase 
loss by 

A0i oss = z/ ma x* - J v(i)dt, (5) 

where v(t) is the derivative of our phase model. We find 
that oscillations typically lose between 1-5 cycles, depend- 
ing on the magnitude of the frequency drift, the duration 
of the oscillations, and the time scale for the frequency 
evolution. One exceptional oscillation, illustrated in Fig- 
ure 6b, drifts by more than 17 cycles with respect to ^ max . 
This oscillation is unusual for its combination of a long 



duration (10 s), a large frequency change (Ais = 4.1 Hz), 
and a slow drift time scale (r = 5.5 s). 

3.7. "Asymptotic" Frequencies of Burst Oscillations 

Aside from the exceptions noted above, the fundamen- 
tal features that led to the adoption of an exponential 
model for the frequency evolution by Strohmayer & Mark- 
wardt (1999) are present in about 70% of burst oscilla- 
tions: the oscillations increase in frequency when they are 
first observed, but appear to reach stable frequencies on 
the time scale of seconds. Several authors have noted that 
the asymptotic frequencies of burst oscillations detected 
several years apart are remarkably similar (Strohmayer 
et al. 1998; Strohmayer & Markwardt 1999; Muno et al. 
2000). Therefore, we measured the maximum frequencies 
that the oscillations reached according to our best fit mod- 
els, in order to determine whether the apparent saturation 
frequencies of the oscillations arc stable. The 1-a uncer- 
tainties on the maximum frequency were found by a search 
through x 2 space (e.g Lampton, Margon, & Bowyer 1976). 
In those instances where the best available fit was statis- 
tically unacceptable, we inflated the uncertainties on the 
phase measurements to for the reduced x 2 to equal 1 be- 
fore the search in x 2 space. The maximum frequency was 
constrained to occur after a minimum, as in several cases 
the absolute highest frequency was observed at the start 
of the oscillation train (e.g Miller 2000, and Figure 3a). 
Our measurements are consistent with those found using 
the Z 2 statistic by Strohmayer et al. (1998), Strohmayer 
& Markwardt (1999), and Giles et al. (2002). We com- 
pare our results for 4U 1636—536 with those of Giles et al. 
(2002) in Section 4.2 

In Figure 9, we plot how the maximum frequencies vary 
with time. To indicate those oscillations where the fre- 
quency has apparently saturated, we have indicated which 
ones required at least an exponential or a 3rd order polyno- 
mial phase model by diamonds. The measurements taken 
from first or second order polynomial models are plot- 
ted only with error bars. We have calculated the frac- 
tional standard deviation in the maximum frequencies, 
ov/(^max), for those sources with more than two values, 
and listed them in Table 3. For Aql X-l and KS 1731-260, 
the frequency evolution appeared to saturate during two 
bursts each, so we report the fractional difference between 
the measurements. The deviations in these frequencies 
ov/ (^max) are quite small, less than one part in 1000. 

4. DISCUSSION 

The frequency evolution of burst oscillations can be de- 
scribed by smooth phase models in all but 30% of the 
pulsations that we have studied (Section 3.5; see also 
Strohmayer 2001). The discrepant cases appear as phase 
jumps in Figure 7. If the oscillations result from brightness 
asymmetries on the surface of the neutron star, it appears 
that some of the following occur: (i) there are multiple 
anisotropics or modes that propagate with slightly differ- 
ent velocities relative to the surface, which results in two 
simultaneous signals that interfere with each other (see 
also Section 3.2), (u) the velocity of the pattern changes 
suddenly on time scales of 0.25 s, or (Hi) the anisotropies 
change their longitude in less than a tenth of a second, 
resulting in sudden phase jumps of about 0.1 cycles. 
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Fig. 9. — The maximum frequencies from Table 1 plotted as the fractional deviation from the mean as a function of time. Measurements 
from oscillation trains which appear to saturate, and are therefore more likely to be indicative of an "asymptotic" frequency, are marked with 
diamonds, while those which do not are simply reported with error bars. 



On time scales of years, the oscillations are stable to a 
few parts in 1000. This is consistent with the Doppler 
shifts that are expected from the binary orbit. Since this 
raises the possibility of measuring a mass function for the 
companion star, we now examine this frequency dispersion 
in detail. 

4.1. The Stability of the Maximum Frequencies 

The change in the apparent spin frequency due to the 
orbital motion of the neutron star is given by 



(^max) pV \ Mx + Me) 2 / 3 

where M x and M c are the mass of the neutron star and 
companion respectively in solar units, Ph r is the binary 
orbital period in hours, and i is the inclination of the or- 
bital plane normal to the line of sight. Assuming a star 
with solar density that over-fills it Roche lobe, the mass 
and orbital period are related by M c w O.llPhr (Frank, 
King, & Rainc 1995). For P hr = 5, sini = 0.87, M x = 1.4, 
and M c — 0.5, the half- amplitude of the Doppler varia- 
tion would be Aj/ max /(j/ max ) = 4 x 10~ 4 . This is a factor 
of 2-5 smaller than the dispersions from 4U 1636—536, 
4U 1702-429, and 4U 1728-34 (Figure 9), but it is of the 
right order of magnitude. 

Only three oscillation sources have known orbital 
ephemerides: MXB 1659-298 (Wachter, Smale, & Bailyn 
2000), Aql X-l (Welsh, Robinson, & Young 2000), and 
4U 1636-536 (Giles et al. 2002). We did not measure 



enough oscillations from MXB 1659-298 and Aql X-l to 
constrain the parameters of the binary orbit. We have 
measured eleven oscillations that appear to saturate from 
4U 1636—536 (the diamonds in Figure 9), so we have at- 
tempted to fit a sinusoid to the maximum frequencies as 
a function of orbital phase (Figure 10). The best-fit si- 
nusoid, allowing the reference phase to vary, has a \ 2 of 
221 for 8 degrees of freedom. Therefore, sinusoidal or- 
bital modulations cannot cause all of the observed dis- 
persion. In Figure 10 we plot a sinusoid with an am- 
plitude equal to the standard deviation of the points fit, 
ov/(^max) = 7.7 x 10~ 4 . The disagreement between the 
data and the sinusoid is primarily due to three oscillations 
with low maximum frequencies. Two of these lasted for 
more than 5 s, one of which is displayed in the right panel 
of Figure 6 (1999 June 10). We will discuss this further in 
the next section. 

There are no other effects that ordinarily change the ap- 
parent spin frequencies of neutron stars (i.e. pulsars) that 
can produce the dispersion in Figures 9 and 10. Accretion 
torques can alternately quicken or slow the rotation of the 
neutron star, and would provide the next largest effect 
(Frank et al. 1995). The torque on a neutron star can be 
estimated by assuming that its magnetic field disrupts the 
disk at some fiducial radius Rm and removes all of its an- 
gular momentum. At Rm, we have Ifl = M(GMRm) 1 ^ 2 ■ 
We take Rm to be the co-rotation radius, where 
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Fig. 10. — The fractional deviation of the maximum frequen- 
cies from Table 1 are plotted as a function of orbital phase for 
4U 1636—536 (see also Figure 6). The dotted line indicates the ex- 
pected orbital modulation using the ephemeris of Giles et al. (2002). 
The dispersion in asymptotic frequencies can not be explained by 
orbital motion alone. 



27w kHz = {GM/R^) 1 / 2 , so that 
v = 7x 10- 5 MReM^ 3 I 



(7) 

where M is the accretion rate in Eddington units, R e is 
the radius of the neutron star in units of 10 km, M x is its 
mass in solar units, I45 is its moment of inertia in units of 
10 45 g cm 2 , and t'kHz is the stellar spin frequency in kHz. 
The values are scaled to those expected for a burster. For 
a spin frequency of 200 Hz, Equation 7 corresponds to 
a Az/ max /(z/ max ) « 6 x 10~ 7 in one year. Clearly, the 
observed dispersion in the maximum frequencies on time 
scales of a month (Figure 9) is much larger than can be 
explained by accretion torques (Equation 7). 

4.2. Frequency Dispersion in 4-U 1636—536 

Giles et al. (2002) recently examined a larger sample of 
oscillations from 4U 1636—536, and found that there is a 
subset of the burst oscillations whose maximum frequen- 
cies are distributed as a Gaussian about 581.64 Hz with a 
width of 0.08 Hz. We cannot independently test this, as 
the uncertainties that we report on the maximum frequen- 
cies are up to a factor of 5 larger than those in Giles et 
al. (2002). This is because the phase-connection technique 
we used to measure the maximum frequencies accounts for 
the frequency evolution of the oscillations, while the dy- 
namic Z 2 technique Giles et al. (2002) used assumes that 
the frequency is constant during each 2 s interval that they 
test. In our analysis, the maximum frequency is usually 
observed at the end of the oscillation train where the data 
can tolerate relatively large variations in the values of the 
polynomial models. 5 In the analysis of Giles et al. (2002), 
if significant frequency evolution occurs during a 2 s inter- 
val, they measure the mean frequency in that interval, but 
do not consider how much larger the frequency could be 
at the end of the interval. As a result the magnitude of 
our uncertainties, 0.05-0.4 Hz, are often larger than the 
width of the Gaussian of Giles et al. (2002). 



If we only include those oscillations that are part of the 
Gaussian in Giles et al. (2002), we find that the data are 
consistent with no modulation. Fitting a sinusoid to the 
data, we derive a fractional amplitude of (5.0±2.1) x 10~ 4 
(x 2 — 6.7 for 5 degrees of freedom). This is equivalent to 
a Doppler shift of 150±63 km s _1 , which is significantly 
larger than the 90% upper limit to the velocity from Giles 
et al. (2002; 55 km s _1 ). The larger uncertainties that 
we assumed result in larger upper limits on the velocity. 
These systematic uncertainties also present an additional 
hurdle in measuring the orbital motion of the neutron star. 



4.3. Models for the Frequency Evolution 

Three models have been proposed to explain the fre- 
quency evolution of the burst oscillations, all of which as- 
sume that anisotropies develop in the surface brightness 
of the neutron star (Strohmayer et al. 1997; Cumming & 
Bildsten 2000; Hcyl 2002; Spitkovsky et al. 2002). The 
simplest model was proposed by Strohmayer et al. (1997), 
who suggested that the oscillations originate from hot re- 
gions on a burning layer that expands and decouples from 
the neutron star at the start of the burst. The oscilla- 
tions are observed when the burning layer begins to con- 
tract, causing them to increase in frequency as the layer 
re-couples to the neutron star. However, recent calcula- 
tions assuming that the burning layer is rigidly rotating 
predict a frequency drift a factor of three smaller than 
that observed (Cumming et al. 2002). Moreover, Cum- 
ming & Bildsten (2000) predicted that there should be 
a fractional dispersion of at most Ai>/i> « 10~ 4 in the 
observed asymptotic frequencies, if the height of burning 
layer is smaller when it re-couples to the star than when it 
decouples. This is a factor of 10 smaller than the residual 
dispersion from 4U 1636—536 in Figure 10. 

As a result, it is important to consider other models for 
the frequency evolution of burst oscillations. Spitkovsky et 
al. (2002) have suggested that a brightness contrast could 
result from a hydrodynamical instability in a gcostrophic 
flow, similar to Jupiter's Great Red Spot, while Hcyl 
(2002) has proposed that Rossby waves propagating retro- 
grade on the neutron star could produce traveling bright- 
ness patterns on the stellar surface. The frequency evolu- 
tion in both models is caused by changes in the velocity 
at which the pattern propagates around the star. The ve- 
locity in turn is sensitive to the temperature, density, and 
vertical structure of the surface layers. Thus, variations 
in the properties of the surface layers near the end of a 
burst can naturally explain the observed dispersion in the 
asymptotic frequencies of the oscillations (Figures 9 and 
10). On short time scales, instabilities in either mecha- 
nism conceivably could cause the apparent phase jumps in 
Figure 7. Moreover, Heyl (2002) also has predicted that 
Rossby waves with different numbers of radial nodes could 
produce simultaneous signals at multiple frequencies, as 
are seen in rare cases (Section 3.2; compare Section 3.5). 
These models are promising, although how these mecha- 
nisms produce brightness variations on the stellar surface 
remains to be studied in detail. 



5 On the other hand, we may be over-estimating our uncertainties when we model the oscillations with higher-order polynomials, which tend 
to have large derivatives at their endpoints. A more physical model might assume that the frequency change in the oscillations is much smaller 
at the endpoints. However, the oscillation that appears 4 Hz from the main signal at the end of the burst on 1999 Apr 14 from MXB 1659—298 
(Wijnands et al. 2001) makes us reluctant to apply this assumption in our analysis, and leads us to assume larger uncertainties. 
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5. CONCLUSIONS 

We have examined 68 burst oscillations observed in 159 
bursts from 8 different sources, and have modeled the fre- 
quency evolution of 59 oscillations from 6 sources (Ta- 
ble 2). The most notable result is that there appears to be 
slight instability in the mechanism generating the oscilla- 
tions on time scales of both seconds and years. About 30% 
of the oscillation trains do not exhibit smooth frequency 
evolution, as both low order polynomials and exponential 
models are statistically inconsistent with the data (Fig- 
ure 7). The possible explanations are that (i) two sig- 
nals may be present simultaneously, invalidating our as- 
sumption that there is a single signal in our analysis, (ii) 
there are discrete phase jumps that occur on time scales 
less than a second, and (in) the frequencies of the oscil- 
lations shift dramatically on time scales of 0.25 s. On 
longer time scales, the maximum frequencies observed dur- 
ing these oscillation trains exhibit a fractional dispersion 
of Azw/(iw) < 4 x 10- 3 (Figure 9). For 5 of the 
sources studied, this dispersion could be consistent with 
Dopplcr shifts of signals originating from rotating neutron 
stars as they orbit the centers of mass of their binaries. 
However, the dispersion in the maximum frequencies in 
4U 1636—536 is uncorrelated with its known orbital pe- 
riod (Figure 10). 

We now are able to accurately characterize the general 
properties of the frequency evolution in our large sample 
of bursts: 

1) The frequency of the oscillations drifts by up to 1.2% 
(Figure 1 and 2), but generally reaches a stable value on 
the time scales of seconds. 

2) Both spin-down in the oscillation train and obvious 
examples of simultaneous signals separated by around 1 
Hz are relatively rare, occurring in 3 and 2 bursts out of 
65 respectively (Figure 3 and 4). 

3) The amounts of the frequency drifts appear larger 
when oscillations are observed earlier in the bursts (Fig- 
ure 8), but are uncorrelated with the durations of the os- 
cillations. This implies that the starts of the bursts set 
in motion the eventual changes in the oscillation frequen- 
cies, as opposed to some mechanism that operates when 



the oscillations themselves appear. 

4) Photospheric radius expansion appears to temporar- 
ily interrupt oscillations. All of the 16 oscillation trains 
that appear continuously from the rise to the tail of the 
burst fail to exhibit radius expansion. When radius ex- 
pansion does occur, oscillations are often seen in the rise 
and/or the tail of the burst, but rarely during the radius 
expansion (Figure 5). 

5) Aside from the above points, the time scale and mag- 
nitude of the frequency evolution does not appear corre- 
lated with other properties of the burst, such as its decay 
time scale, peak flux, or fluence. 

6) In about 60% of the oscillation trains, the frequency 
evolution displays inflections that cannot be described 
by an exponential model. Two illustrative examples are 
shown in Figure 6. Thus, in addition to the upward fre- 
quency drift, the frequency tends to wander by on order 
0.1 Hz in a random manner. 

Several models have been proposed in which the oscilla- 
tions originate from an anisotropy in the emission from the 
neutron star's surface. Each models the frequency evolu- 
tion as originating from different mechanisms (Strohmayer 
et al. 1997; Cumming et al. 2002; Heyl 2002; Spitkovsky 
et al. 2002). The apparent lack of stability in the underly- 
ing clock producing these oscillations favors the models in 
which the brightness asymmetry originates from hydrody- 
namic instabilities (Spitkovsky et al. 2002) or modes ex- 
cited in the neutron star ocean (Heyl 2002). Future stud- 
ies of the amplitude, harmonic structure, energy spectrum, 
and phase lags of the oscillations will further constrain the 
location of the hot spot and probe the atmosphere of the 
neutron star. 

We thank Lars Bildstcn and Fred Lamb for useful 
comments, and the referee for careful reading of this 
manuscript. We also thank Pavlin Savov for developing 
the burst detection software, and Derek Fox for develop- 
ing much of the phase connection software. This work 
was supported by NASA, under contract NAS 5-30612 and 
grant NAG 5-9184. 
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Table 1 

Summary of Burst Oscillations 





J^burst 




Number of 




Source 


(Hz) 


Bursts Osc. Models a 


4U 1608-522 


620 


6 


1 





4U 1636-536 


581 


24 


17 


17 


MXB 1659-298 


567 


15 


5 


3 


Aql X-l 


549 


17 


3 


3 


KS 1731-260 


524 


13 


5 


4 


4U 1728-34 


363 


66 


27 


24 


4U 1702-429 


329 


11 


8 


8 


4U 1916-053 


270 


7 


1 






a The number of oscillations modeled using the 
phase connection technique. 
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Table 2 

Models for the Phase Evolution of Burst Oscillations 



Best 


Number 


Model 


of Oscillations 


1st Order 


9 


2nd Order 


13 


3rd Order 


18 


4th Order 


3 


5th Order 


10 


Exponential 


6 



Table 3 

Dispersion in Maximal Frequencies of Burst Oscillations 

Source av/(f ma x) 

4U 1636-536 7.7xlCT 4 

4U 1702-429 7.5xl0~ 4 

4U 1728-34 9.2xl0~ 4 

KS 1731-260 a 7.3xl0~ 5 

AqlX-l a 1.2xl0~ 4 



Note. — o y is defined as 
the standard deviation of the 
maximum frequencies mea- 
sured in oscillations that ap- 
pear to saturate, unless oth- 
erwise indicated. 

a Only two maximum fre- 
quencies were measured, so 
cr„ is the difference between 
the values. 



